Initialisation
## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin Phi[2] = valeur du noeud Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,90,5),
omega2 = c(0.5,0.1,0.01) )
#=======================================#
t <- seq(60,120, length.out = 10) #value of times
# --- longitudinal data --- #
set.seed(18031997)
data <- NLME_data(G = 20, ng = 8, time = t, fct = m, param = param)
getDim(data)
## G ng N n F.
## 20 8 160 1600 3
Y <- data$obs

SAEM avec simulation par MCMC
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
Phi <- fct_vector(function(sigma2, rho2, mu, omega2)
c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), mu), dim = c(1,1,F.,F.) )$eval
S <- fct_vector(function(eta, phi) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
function(eta, ...) mean(eta^2),
function(phi, ...) apply(phi^2, 2, mean),
function(phi, ...) apply(phi , 2, mean), dim = c(1,1,F.,F.) )
Metropolis Hastings
loglik.phi <- function(x, eta, Phi)
{
id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
sum(Phi%a%1*S$eval(eta = eta, phi = setoffset(x, Phi%a%4), i = 1) + Phi%a%3 * S$eval(phi = x, i = 3) )
}
loglik.eta <- function(x, phi, Phi)
{
id <- c(1,2)
sum( Phi%a%id * S$eval(eta = x, phi = phi, i = id) )
}
SAEM
Initialisation
# --- Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.2 ; para$omega2 <- rep(.1,3)
# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( matrix(rep(0, G*ng), ncol = 1) ),
phi = list( matrix(rep(para$mu, G), nrow = F.) %>% t ) )
Etape simulation et maximisation du SEAM
sim <- function(niter, h, Phih, eta, phi, verbatim = F)
{
M <- length(phi)
eta <- 1:M %>% lapply( function(i)
MH_High_Dim_future(niter, eta[[i]], sd.eta , loglik.eta, phi[[i]], Phih,
cores = 1, verbatim = verbatim))
phi <- 1:M %>% lapply( function(i)
MH_Gibbs_Sampler_future(niter, setoffset(phi[[i]],-Phih%a%4), sd.phi, loglik.phi, eta[[i]], Phih,
cores = 1, verbatim = verbatim )) %>%
lapply(setoffset, Phih%a%4)
list(eta = eta, phi = phi)
}
maxi <- function(S)
{
list(sigma2 = S%a%1,
rho2 = S%a%2,
mu = S%a%4,
omega2 = S%a%3 - (S%a%4)^2 )
}
|
|
sigma2
|
rho2
|
mu1
|
mu2
|
mu3
|
omega21
|
omega22
|
omega23
|
|
Oracle
|
0.0502965
|
0.1040354
|
4.841866
|
89.83036
|
5.020998
|
0.2875514
|
0.1163807
|
0.0098201
|
|
Initialisation
|
0.0545425
|
0.2000000
|
5.763583
|
103.74449
|
5.763583
|
0.1000000
|
0.1000000
|
0.1000000
|
niter <- 100
MH.iter <- 10
sd.eta <- .3
sd.phi <- c(.05, .3, .05)
gg <- simulation_test(sim, Phi, param, 100, Z) %>% lapply(function(z)plot(z[[1]], nrow = 2))


res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi,
burnin = 150, coef.burnin = 2/3, verbatim = 2)
saveRDS(res, 'saem.rds')
[1] “SAEM execution time = 01min 29sec”
Result of the SAEM-MCMC
|
|
sigma2
|
rho2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
|
Real value
|
0.0503
|
0.104
|
4.8419
|
89.8304
|
5.0210
|
0.2876
|
0.1164
|
0.0098
|
|
Estimated value
|
0.0499
|
0.095
|
4.8736
|
89.8251
|
5.0863
|
0.3167
|
0.1535
|
0.0182
|
|
Rrmse
|
0.0079
|
0.087
|
0.0066
|
0.0001
|
0.0130
|
0.1014
|
0.3185
|
0.8520
|
## [1] "SAEM execution time = 01min 29sec"

## $plot_parameter

##
## $plot_MCMC

##
## $plot_acceptation
## TableGrob (2 x 1) "arrange": 2 grobs
## z cells name grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]
